Modeling the elastic deformation of polymer crusts formed by sessile droplet 
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Evaporating droplets of polymer or colloid solution may produce a glassy crust at the liquid- 
vapour interface, which subsequently deforms as an elastic shell. For sessile droplets, the known 
radial outward flow of solvent is expected to generate crusts that are thicker near the pinned contact 
line than the apex. Here we investigate, by non-linear quasi-static simulation and scaling analysis, 
the deformation mode and stability properties of elastic caps with a non-uniform thickness profile. 
By suitably scaling the mean thickness and the contact angle between crust and substrate, we find 
data collapse onto a master curve for both buckling pressure and deformation mode, thus allowing 
us to predict when the deformed shape is a dimple, mexican hat, and so on. This master curve is 
parameterised by a dimensionless measure of the non-uniformity of the shell. We also speculate on 
how overlapping time-scales for gelation and deformation may alter our findings. 

PACS numbers: 46.32,+x, 81.15.-z, 46.70.De 
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I. INTRODUCTION 

The emergence of cost-effective microfiuidic devices al- 
lowing manipulation and control of fluids on microme- 
ter scales has promised a significant new paradigm for 
the manufacturing industry P, • Known as bottom-up 
processing, the required material can be transported to 
its destination in solvent form and deposited onto the 
target substrate using established technologies such as 
ink-jet head fluid ejection. One constraint of this means 
of material transport is that the deposited solution will 
evaporate, inducing non-trivial solvent flow that may dis- 
tribute the solute in an irregular and undesirable manner. 
Understanding the physical mechanisms underlying such 
phenomena are thus crucial to obtaining the desired de- 
posit profile. 

Of concern here is the formation of a glassy 'crust' 
that can appear on the liquid-vapour interface of evap- 
oration polymer [1 H H |l or colloid 0, H solu- 
tions. Such crusts are believed to produce roughened 
surface profiles after drying of thin films ^3 • F° r sessile 
(i.e. pinned contact line) droplets on a wetting substrate, 
the effect is compounded by the outward radial flow of 
solvent during evaporation which, for low-concentration 
solutions, leads to enhanced solute deposition near the 
pinned line 0, IH El 13 H. For the higher concen- 
trations where crust formation occurs, this outward flow 
suggests that the thickness of the crust will vary non- 
uniformly over the surface of the droplet, being thicker 
near the outer contact line and thinner near the centre 
or apex. 

Being elastic, such a crust will deform under the in- 
ternal osmotic pressure generated by continued solvent 
evaporation through its porous surface, and thus the final 
deposited polymer profile will at least partly depend on 
the mechanical properties of the crust. Assuming it to be 
everywhere thin, the relevant area of elasticity required is 
shell theory, an established field with known results for 



many different geometries, including the spherical caps 
relevant here (the droplet surface is spherical, so the ini- 
tial crust profile should be approximately so); see e.g. 
Shilkrut |l5j for a recent survey. However, most of the 
literature is for shells of uniform thickne ss |l6l UtI ITsl ITsj 
or specific, engineered non-uniformities 20]; those with a 
thickness profile similar to that expected for evaporating 
droplets have not been investigated. Furthermore, the 
rich array of deformed shell profiles observed in experi- 
ments UyHSEiH nas n °t, to the best of our knowledge, 
been systematically quantified as a function of the shell 
parameters, even for uniform shells. 

In this paper we describe numerical and theoretical in- 
vestigations into the deformation of closed, elastic spher- 
ical caps with a thickness profile expected of crust for- 
mation during droplet evaporation, namely thin near the 
apex and thicker near the contact line. We restrict our- 
selves here to axisymmetric deformations that preserve 
the axis of symmetry of the shell; asymmetric deforma- 
tions, as sometimes seen in experiments will be the 
subject of a future study. An overview of the problem 
is given in Fig. ^a). Here, two schematic equilibrium 
curves of inward pressure P and change in droplet volume 
AU are given, one monotonic and one non-monotonic 
(here and below we define P and AU to be positive for the 
deformations of interest). The S-shaped non-monotonic 
curve exhibits a buckling instability at a value P c , when 
the shell jumps to an approximately inverted shape with 
a boundary layer. This is known as snap-through, or 
simply snap buckling, and can also be realised by e.g. 
applying a localised load ^3 or long-range attractive 
force |2l|. If the pressure is subsequently decreased, a 
lower critical pressure is reached when the shell jumps to 
the original solution branch. This instability is known as 
snap-back buckling. The two differing critical pressures 
gives a hysteresis curve, the integrated area of which cor- 
responds to the combined energy dissipated during both 
dynamic buckling events as the shell is damped to the 
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n=1:dimple 



n=2: mexican hat 



(b) 



n=3: crater 



FIG. 1: (a) Schematic equilibrium P — AV curves, one that 
buckles (with a dashed unstable regime), showing a hystere- 
sis curve for cyclically-varying P, and one that does not. 
Sketches of possible shell shapes are also shown, (b) Lat- 
eral cross sections of pre-buckled shells (black) compared to 
the undeformed shape (grey), showing the first 3 deformation 
modes. For the crater mode, the turning points of the profile 
are indicated by arrows. 



static, equilibrium curve. 

The above discussion assumes the pressure P is being 
controlled. In fact, for the evaporation problem it is more 
natural to consider a slowly varying AV , which may ex- 
hibit different limit points; indeed, the simple schematic 
curves given in the figure are always stable under con- 
trolled volume. Nonetheless we shall here control P, and 
focus much of our attention on the snap-through buck- 
ling point. The principle reason for doing this is that we 
are most interested in the deformation of the shell once 
non-linear effects have first become established, and the 
snap buckling point is a convenient point to focus on. 
Indeed, as shown below, the mode of the shell remains 
essentially unaltered from when non-linear effects first 
become significant to the buckling point itself, so we are 
effectively probing a broad range of the P — AV curve. 
Furthermore, controlling P permits a closer correspon- 
dence with the spherical cap literature, for which a con- 
trolled volume is not typically considered. It also avoids 
the need for more sophisticated (and time-consuming) 
incremental-iterative algorithms p^ . 

This paper is arranged as follows. In section [H] the 
model and simulation method are described. Results are 
then presented in section ITTT1 where it is shown that data 
for P c as a function of the undeformed shell geometry 
can be made to collapse onto a single curve when scaled 
accordingly. The same scaled parameters also determine 
the shape of the shell, as quantified by the 'mode number' 
n (see Fig. ^b) for some examples) . By parameterising 



shells of non-uniform thickness in a convenient way, we 
show that the effect of non-uniformity is to distort these 
master curves, without altering the gross underlying be- 
haviour. This is supported by the scaling analysis of 
section ITVl Finally, we return to the droplet evaporation 
problem in light of our findings in section IVI 



II. MODEL AND SIMULATION METHOD 

To focus on the behaviour of elastic shells of non- 
uniform thickness, we shall here reduce the complex pro- 
cess of crust formation into a simple, idealised form. 
In particular, modulation of the liquid- vapour interface 
during crust formation, crust thickening after its initial 
genesis, and the possibility of temporally-evolving intrin- 
sic curvature are all ignored. In essence we assume that 
the all of the polymer rapidly gels at the surface of the 
droplet, and is initially in a stress-free state. Since liq- 
uid droplets of size less than a few mm have spherical 
surfaces |23| . our initial, undeformed shell configuration 
is taken to be a spherical cap with radius of curvature R 
and a contact angle (the angle the shell makes with the 
substrate) 9q. This curved undeformed state identifies 
the crust as a shell, distinct from a plate which is flat 
when not deformed |24j . 

The static configuration of a thin shell can be purely 
defined in terms of the geometry of the mid-plane surface 
S and a (possibly variable) thickness h, defined such that 
the outer and inner surfaces lie at distances ±/i/2 normal 
to S. As we are only concerned here with quasi-static 
deformations, it suffices to find shell configurations that 
minimise U — PAV, with U the elastic strain energy and 
— PAV is the work done by the pressure P (in our sign 
notation). It is convenient to partition IA into two parts, 
a stretch term 
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where Uij is the strain of the mid-plane, in spherical 
coordinates in which 9 G [0, 9q] is the azimuthal and 
4> £ [0, 27t] the polar angles, resp., E is the Young's mod- 
ulus of the material, and v its Poisson ratio. In addition, 
there is a bending energy term 
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where Ky the change in curvature of the shell meridian 
from the undeformed value l/R p5ll26j|. Terms propor- 
tional to ug§ or Kg^ vanish for axisymmetric deformations 
and have been dropped. 

The expressions Q and J3J) allow for E, v and the 
thickness h to vary over the surface. We shall here only 
consider variations in thickness that respect the axis of 
symmetry of the shell, i.e. h = h(9). For simplicity, a 
two-parameter quadratic form for h(9) is employed, 



3 



h{9)=A({h),p)+B({h),p)9 2 . (3) 

Here, the constants A and B are completely specified 
by the mean thickness (h) and a dimensionlcss non- 
uniformity parameter p — [h(9o) — h(0)]/h(9o). In this 
way, p = corresponds to uniform shells, and p > to 
shells that are thinner near the apex than the contact 
line. The apex thickness vanishes in the limit p —* 1 _ . 

Equilibrium configurations S for given shell parame- 
ters and pressure were found by numerically minimis- 
ing the total energy using the BFGS variable metric 
method This was performed on the displacement 

vector fields u r (9) (radial displacement) and ug(9) (an- 
gular displacement), where the 9 range was discretised 
into 50 — 100 segments as the situation dictated. The 
midplane strain tensor mj was found by evaluating the 
full, (geometrically) non-linear expression 2ua = diUj + 
djUi + (diUk)(djUk) in spherical coordinate s 1241 . For sim- 
plicity, only the linear expressions for «y |26j were used, 
since intuitively it is non-linearity in the in-plane strains 
that controls buckling. The severity of this simplification 
was checked by comparing simulation data to known re- 
sults for spherical caps (i.e. the buckling pressure P c for 
uniform h p"?l Il9jh for which the agreement was good. 

The results presented below are for a Poisson ratio 
v = 0.3 and clamped boundary conditions, in which the 
angle the shell makes with the contact line is fixed. For 
robustness, we have also simulated a representative sam- 
ple of systems with v = 0.5 (incompressible), v = and 
v = —0.5, and found no qualitative change in the results, 
merely a modulation of the master curves (discussed be- 
low) amounting to no more than 10%, although inspec- 
tion of the shell equations IIIL'I) suggests values of v close 
to —1 might produce more significant deviations. For 
hinged (i.e. no angle constraint) boundary conditions 
the effect was larger, reducing A cusp (again, see below) 
by « 40%, and similar sized changes in the critical pres- 
sure. Furthermore the master curves for non-uniform 
shells were primarily shifted in pressure, only weakly in A. 
Considerations of the microscopic nature of the crust- 
substrate interaction suggest clamped boundaries to be 
more realistic, and so we focus our attention on these. 

III. RESULTS 

For small pressures P, the shell displacements are lin- 
ear in P and the elastic strain energy is simply propor- 
tional to P 2 . Furthermore the shape of the shell always 
takes the lowest mode consistent with the boundary con- 
ditions, namely the n — 1 'dimple' form. Higher modes 
only arise when the non-linear terms in the membrane 
strain become important. This is demonstrated in Fig. El 
which shows the radial displacement as a function of 9 
for different pressures P from up to the buckling point. 
Defining the mode n to be the number of turning points in 
the radial displacement u r (9) in the range < 9 < 9q, it 
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FIG. 2: Example of the evolution of the shell deformation 
with pressure. Each curve corresponds to the outward ra- 
dial displacement u r (8), scaled by the uniform shell thickness 
(h) = h = 5 x 10" 3 i?, for different pressures P. The final 
value of P is just prior to buckling. In this example, 9q = 0.5, 
v = 0.3 and E is the Young's modulus of the shell material. 



is evident that n = 1 for small P but takes a higher value 
just before buckling, n = 4 in this example. Note that 
n does not pass through any intermediate states n = 2 
or n — 3; the mode n of the shell deformation is fixed 
once non-linear effects become established, only the am- 
plitudes continue to vary with P. This is typical of the 
behaviour observed in the parameter range studied. Thus 
when below we present results for n just prior to buck- 
ling, it should be understood that this same mode exists 
for all prior P down to the linear regime. 

The physical mechanism behind the snap-buckling is 
as follows. In the limit /&—►(), the bending energy 
i^bend © becomes arbitrarily small compared to the 
stretching energy £4tretch (QJ and so can be ignored (un- 
less the shell deforms inextensibly, which is not possible 
for this geometry |2a l28|). Now observe that if the de- 
formed shell takes the shape of an inverted spherical sur- 
face with the same radius of curvature R, the midplane 
strains Uij, and hence i4tretch> vanish. In between this 
inverted shape and the undeformed non-inverted shape 
lay flatter configurations with a high £4trctch, which col- 
lectively play the role of an energy barrier between the 
two preferred states. As P is varied, this energy land- 
scape 'tilts' until the system hops between the two local 
minima, i.e. it snap-buckles. 

In contrast to the stretching energy, the bending en- 
ergy Ubend increases monotonically, even when the shell 
becomes inverted. Hence as (h) increases and £4>cnd be- 
comes more significant, the depth of the energy well cor- 
responding to the inverted shape becomes smaller, until 
at some point it vanishes. When this happens, there is 
no buckling and the shell continuously and smoothly de- 
forms for slowly varying P (i.e. the P — AV curve is 
monotonic). Thus it should be possible to identify some 
dimensionless combination of parameters, including (h), 
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FIG. 3: Plot of shell shape just prior to buckling (as given in 
the key) as a function of contact angle 9q and shell thickness 
normalised to the radius of curvature, h/R, for a uniform 
h{9) = (h). The 'blockiness' is due to a large data interval for 
#o of w 41%. The solid circles denote the thickness at which 
buckling is first observed, and the dashed line has the slope 
1/2 as shown. 



that determines whether or not the shell buckles; for (h) 
too large, there should be no instability. This insight is 
confirmed the simulation results, as we now describe. 



A. Critical pressure P c and mode n 

The variation of the mode n prior to buckling with 
both contact angle 6*o and uniform shell thickness h/R is 
given in Fig. EJ As expected, for each 9q there is a critical 
thickness that separates shells that buckle (thin shells) 
from those that don't (thick shells). For small angles 
this line is well approximated by (h)/R oc #g, which can 
be attributed to a crossover from a stretch-energy domi- 
nated regime to a bending-energy dominate one; see for 
instance the scaling theory of [2J] or section llVl Further- 
more, the mode n increases the further one moves into the 
buckling regime, and again the crossover between modes 
seem to lie on lines (h)/R oc 9q, an observation that is 
quantified below. Thus dimpling, mexican hats and so on 
are just the beginnings of a series of modes that increases 
apparently without bound. 

In the so-called shallow shell regime 6*o 1, it is cus- 
tomary to describe the shell geometry by the dimension- 
less quantity 



A =[12(1- f)\ 



(h) 



(4) 



Intuitively, small A corresponds to thick and/or flat 
shells, which behave something like flat plates, and large 
A to thin and/or steep shells which deform as curved 
shells. In addition, the pressure is normalised to the 
buckling pressure for a full sphere psj , 
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FIG. 4: Buckling pressure P c scaled to the full sphere solution 
psphere ( U pp er p ane i) anc j the mode just prior to buckling 
(lower panel) versus A for shells of uniform thickness h and 
contact angles Oq given in the key. The data range for each 
6<o differ, but all extend down to A cusp . 
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When plotted in terms of these dimensionless variables, 
both P c and n collapse onto a single curve, as shown in 
Fig. Although the highest 9q were omitted from this 
plot, deviations only amount to around « 15% for 9 = 

0. 8 rads « 45°, suggesting the so-called 'shallow shell' 
limit provides a decent approximation for even somewhat 
steep contact angles. 

It is apparent from this figure that q c f=a 1 for A > 1, 

1. e. the buckling pressure approaches the full sphere so- 
lution for large values A. Furthermore, the mode n ~ A 
to good approximation, which is also confirmed by the 
scaling argument given later. Plotting both sets of data 
on the same axes as here also shows a correlation be- 
tween the features of the curves, namely that the 'kinks' 
in the pressure data correspond to increases in the mode. 
Thus deviations from the complete sphere solution de- 
pend on the mode of deformation. As a final observation, 
note that a minimum A cusp » 3.2 for buckling to occur 
was observed for each data curve (the label 'cusp' is ex- 
plained below). This, and the shape of the P C (A) curve, 
lie within 5% of known results from more rigourous shell 
simulations 15] , justifying the simplifying choice of linear 
bending expressions in the simulations. 

Introducing a variable shell thickness through the non- 
uniformity parameter p > does not alter the data col- 
lapse described above, but rather shifts and distorts the 
master curve. P c versus A for different values of p is given 
in Fig. [SI in terms of the same dimensionless quantities 
as before. There is again good data collapse, but now 
onto distinct curves, with one such curve for each p. As 
p increases, the curves move to lower pressures; in words, 
shells that are thinner near the apex than the contact line 
are weaker than uniform shells with the same mean thick- 
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?i = (12(1-v 2 )) 1/4 (e 2 R/<h» 1/2 

FIG. 5: Buckling pressure P c scaled by the full sphere solu- 
tion, for different values of the non-uniformity parameter p 
given in the key. For each value of p, data for 80 = 0.1, 0.2 
and 0.4 has been plotted. To guide the eye, an arrow have 
been added to each line to indicate the kink corresponding to 
the n = 1 «-> n — 2 mode transition. 

ness. Furthermore they move to higher A, so that A cusp 
increases and a broader range of shells will not buckle 
(i.e. have A < A cusp ). Note that although the curves 
for different p share the same essential features, includ- 
ing kinks when the mode number changes, they cannot 
be collapsed onto a single 'super-master' curve by simple 
scaling. Note also that the relationship n ~ A still holds 
(data not shown). 

An alternative way to probe non-uniform shell thick- 
ness is to hold A fixed and instead vary p. Since the mas- 
ter curves tend to drift to lower q and higher A with p, in- 
creasing p should act to decrease the mode number (pos- 
sibly removing buckling) and lower the critical pressure. 
This is precisely what is observed; as evident in Fig. 
P c is a decreasing function of p, and furthermore mode 
reduction or loss of buckling can occur at finite values of 
p > 0. Thus a non-uniform shell, rather than increasing 
the mode number (as perhaps might be expected), on 
the contrary acts to decrease the mode. Quantifying the 
dependency on p is difficult as it depends on the details 
of the family of master curves parameterised by p, but 
we note that the data for the 9q = 0.4 case in the figure 
appears to scale as P c (p) ~ (1 — p) 2 for small 1 — p. We 
have no explanation for this at present. 

B. The onset of buckling: Cusp scaling 

Although evaporation is clearly a unidirectional pro- 
cess, it is nonetheless instructive to consider cyclic vari- 
ations around the buckling point since, as evident from 
Fig. ^ such a protocol probes the hysteresis curve and 
allows greater insight to be gained into the nature of 
buckling singularity. A sample of hysteresis curves for a 




FIG. 6: P c as a function of the non-uniformity parameter p, 
scaled to the uniform thickness p — case, for mean thickness 
and contact angle given in the key. The 80 — 0.2 shells cease 
to buckle for p > 0.45. Also, the radial displacement u r (6) 
just prior to buckling for the two shells either side of the 'kink' 
in the 80 = 0.4 line are given, showing the decrease in mode. 
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FIG. 7: Strain energy versus pressure P for 80 = 0.2 and the 
A given in the key, in units such that R — E — 1 . The vertical 
arrows denote critical pressures for snap-through buckling un- 
der increasing P (upwards-pointing arrows), and snap-back 
buckling when P is subsequently decreased (down arrows). 
(Inset) Difference in the two critical pressures AP C and the 
deflection of the apex just prior to buckling Au,-(0) as a func- 
tion of e = A - A cusp with A cusp w 3.231. The solid lines have 
slopes of 1/2 and 3/2 as shown. 



range of A near A cusp are given in Fig. The 'width' 
of the hysteresis loop was taken to be the difference be- 
tween the two critical pressures for snap-through and 
snap-back buckling, AP C . For the 'height', either the dif- 
ference in strain energy AU, or the difference in apex dis- 
placement Au r (9 = 0), between pre-buckled and post- 
buckled states can be used; both quantities give the 
same scaling. As evident from the figure, the hystere- 
sis loop vanishes continuously as A — > A cusp , according to 
AP C - (A - A cusp ) 3 / 2 and Au r - (A - A™ 8 ?) 1 / 2 . 
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These exponents are characteristic of a cusp singular- 
ity [2^, |3(| , and indeed the scaling theory of section II VI 
which explicitly has a potential energy function capable 
of generating a cusp, also predicts the same exponents 
(see below for details) . This explains our use of the su- 
perscript 'cusp ' when referring to this singularity. The 
picture is as follows: for A > A cusp , the system encoun- 
ters two buckling transitions on suitably varying P, one 
corresponding to forward {'snap-through') buckling, the 
other to reverse or 'snap-back ' buckling. As A cusp is ap- 
proached, these two transitions approach until 'annihilat- 
ing' at a higher-order singularity, namely the cusp point. 
For A < A cusp there is no buckling and the shell deforms 
continuously for all pressures. 

IV. SCALING THEORY 

It is possible to construct a scaling argument to de- 
scribe and compare to the simulations results. This is 
similar in style to that in Landau and Lifshitz |24|, but 
with the addition of explicit non-linear terms for the 
stretching energy, and a new variable to mimic the mode. 
We shall also present some simple, linear calculations for 
non- uniform shells. 



A. Variation with mean thickness: Buckling 

The approximations employed in setting up the scal- 
ing theory are drastic but far-reaching. Firstly, only the 
radial component of the displacement u r is considered; 
the angular component ug, which corresponds to in-shell 
motion and is anyway small, is simply neglected. Sec- 
ondly, the angular dependence u r {&) is subsumed into 
a characteristic inward displacement C ~ -ti r , with an 
unknown (but positive as defined here) dimcnsionlcss 
prefactor. We further define a dimcnsionless variable 10 
through the characteristic derivatives du r /d6 ~ —(u/9o 
and d 2 u r /d0 2 ~ —((uj/9o) 2 . We shall below identify 10 
with the mode number n, although u is a continuous vari- 
able and as such potentially contains more information 
than the discrete n. 

In terms of the dimcnsionless quantities A already de- 
fined in (0J, p = P/(E9q) and x — (,/{h), the characteris- 
tic energies per unit surface area (including the work W 
done by the pressure) can be derived from lllll'l) . 

Stretch = E{h)9llCl—- £3—^- + Cl^g- > 

A7J — E(h)6QX 2 , . 2 . 3 , 4 , 

OMbend — \'^2 LU + VLy^UJ + C 4 W j- 

SW =-C w pEe%(h)x (6) 

where all of the C;'s are unknown, but positive, dimen- 
sionless constants. As in the simulations, the non-linear 
strain un but linear curvatures have been used. These 



equations must be supplemented by the constraint that to 
is bounded below by some positive 0(1) constant, which 
is required for consistency with the fixed boundary con- 
ditions. The equilibrium x and lu for a given p can be 
found by minimising the total energy. 

Consider first the linear response x <C 1, for which 
terms in x 3 and higher can be ignored and the energy 
function is quadratic in x. Then (jSJ predicts a negative 
ui, which is not consistent with the boundary conditions, 
so we fix a; to be an 0(1) constant and minimise for x 
alone. In the bending-dominated regime A <C 1, when 
<K/4tretch can De neglected, x ~ pX 8 or £ ~ POqR 4 /E{h) 3 . 
Conversely, in the stretching-dominated regime A > 1, 
x ~ pX 4 or £ ~ PR 2 / E(h). The crossover between these 
two regimes lies at A = 0(1), i.e. 6 ~ (h)/R, which co- 
incides with the onset of buckling observed in the simu- 
lations. This confirms the intuitive picture that buckling 
arises when the stretching term dominates. 

Let us now turn to non-linear effects and buckling. For 
x = 0(1), the energy © is quartic and thus admits 1 or 
3 equilibrium solutions. This is most clearly seen in the 
stretching regime A> 1, when the strain energy is dom- 
inated by ££4tretch an d all 3 roots exists, two minima at 
x = and x ~ X 2 ju) 2 and an unstable local maximum in- 
between. As p increases, the root initially at x — moves 
to higher x until it and the local maximum annihilate 
each other; this is the critical pressure when the system 
jumps to the second minimum. It is straightforward to 
show that, just before this buckling event, x c — 0(1), 
lo c ~ A and p c ~ A~ 4 , or in non-normalised variables, 
Cc ~ (h), uj c ~ 9 ^R/{h) and P c - E{h) 2 /R 2 . This is in 
agreement with the simulation results discussed above if 
we regard w as a continuous analogue of the mode num- 
ber n. It should be mentioned here that the second root 
and local maximum do not exist for all values of the Cf. 
Since it is beyond this scaling theory to determine the 
Cf , we must simply assume they take values such that 3 
roots exists for A ^> 1 . 

For A = 0(1), the size of the jump to the second min- 
imum at p c is smaller than for large A, and vanishes at 
a finite value A cusp . The corresponding critical pressure, 
p£ usp , and A cusp describe a point in parameter space for 
which all three equilibria (including the unstable maxi- 
mum) coincide. To analyse scaling behaviour around this 
point, we define e = A — A cusp and again assume that u) 
is a fixed, 0(1) constant (this is consistent with the sim- 
ulations, for which the lowest mode solution was always 
observed around the critical A) . By expanding the energy 
around the cusp point in e and preserving only the low- 
est order terms, it is straightforward to derive the scaling 
observed in the simulations, namely that the difference 
between the two critical pressure (for snap-through and 
snap-back buckling) scales as Apf ~ e 3 / 2 , and the cor- 
responding configurations scale as Axf ~ e 1 ! 2 . These 
exponents and the energy function 10 describe a cusp 
catastrophe [2^.l3p| . as suggested by our notation. 
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B. Non-uniform thickness 

The nature of scaling theory means that it is not 
well suited to describing spatially- varying quantities. 
Nonetheless a simple argument, in which bending energy 
and non-linearities are ignored, can be treated analyti- 
cally, as we now describe. For p > 0, the apex is thinner 
than the 'wing' region near the contact line, and might 
be expected to deform more. We therefore define two 
characteristic deformations, Ca and £ w , for the apex and 
wing regions respectively. These two regimes are consid- 
ered to be sharply separated at an angle 8\ < 8 , which 
defines the third variable of this model. For any given 
pressure P, the values of Ca, Cw and 8\ are found by min- 
imising the linear strain energy restricted to stretching 
terms, 



U = C a E(l / d6h(8)sm8 + C h E(l / d8h(8)sm8 

JO JBx 

~C ap P( a R 2 [ d8 sind - C wp P( w R 2 [ ° A8 sin 8 

JO JBi 

(7) 

where as before the Cj's are 0(1) dimensionless con- 
stants. Substituting the explicit quadratic form for 
h(8) © is most easily treated by considered small p and 
p w 1 separately, and again we assume shallow shells 

<1. 

For p <C 1, i.e. almost-uniform shells, 8\/8q « 
(1 — p/8)/\/2 and that the deformation around the apex 
is enhanced with respect to the wing by a small amount 
Ca/Cw ~ 1 + p/2, to leading order in p. The criti- 
cal pressure, which can be estimated in linear theory 
as when the shell becomes flat, is similarly reduced, 
P c (p)/P c (0) w (1 - p/4). Although this correctly pre- 
dicts an initially linear decrease in P c with p, it is clear 
from the simulations that the slope depends on the shell 
parameters rather than taking the fixed value —1/4 as 
suggested here. 

Repeating the procedure for p w 1 and expanding in 

1 — p -C 1 reveals a more drastic modification from the 
uniform case. Now 8i/8 a (1 — p) 1 / 4 and Ca/Cw ~ 
(1 — p)~ x / 2 , so a small region around the apex becomes 
strongly deformed. The small exponent of 1/4, however, 
suggests p must be very close to 1 for this regime to 
be clearly observed. The critical pressure vanishes as 
P c (p)/P c (0) w (1 — p) a with a = 1, whereas the simula- 
tions instead suggest a — 2 as already described, again 
demonstrating this argument is at best capable of pro- 
ducing qualitatively correct predictions only. It is, how- 
ever, not clear at this time how to improve the theory 
with regards non-uniform shells while retaining the ap- 
pealing clarity of the scaling approach. 



V. DISCUSSION 

The findings outlined in this article demonstrate that 
the effects of non-uniform shell thickness on the defor- 
mation and buckling of elastic caps can be handled in 
a systematic manner, once the mean shell thickness (h) 
has been factored out of the thickness profile h(8). For 
the quadratic profile considered here, the dimensionless 
non-uniformity parameter p was shown to determine the 
behaviour of the system, in that each p corresponded 
to a different master curve once the various quantities 
have been normalised in a suitable manner. Based on 
this finding, we hypothesise that any non-uniform thick- 
ness profile h(8) = (h)f(8;8o) will produce similar be- 
haviour, with master curves parameterised by f(8;&o). 
The overall form of the curve, such as kinks correspond- 
ing to changes in the mode, is also expected to remain 
the same, with only quantitative details differing. 

The issue then becomes one of determining the actual 
thickness profile generated by droplet evaporation of ses- 
sile polymer or colloid solutions. This is far from triv- 
ial; crust formation involves a number of physical pro- 
cesses [HI . including vapour diffusion, solvent flow and 
gelation, any and all of which may overlap in time with 
the elastic deformation considered in isolation here. Until 
such a time that these effects are properly quantified for 
the droplet geometry, it is difficult to provide any definite 
predictions. Nonetheless some broad statements can be 
made based on our findings. For instance, we expect that 
shells that are thinner near the apex than the base will 
generally have a lower mode number n that uniform shells 
of the same mean thickness. This, perhaps ironically, 
suggests that, if a uniform deposit profile is required, it 
may be advantageous to have a non-uniform thickness, 
since lower mode numbers n correspond to flatter final 
shapes. Also such shells are weaker, in that the critical 
pressure is lower, suggesting they will reach non-linear 
region of the P-AV curve comparatively quickly. 

Notwithstanding the difficulties associated with the 
full evaporating droplet problem, it may be worthwhile 
to briefly speculate here on the effects of ongoing crust 
formation during deformation. Supposing that polymer 
is added to the crust in an initially unstressed state, the 
effects of continued gelation will be to reduce the mean 
shell stress, and alter the intrinsic curvature towards the 
current value. Given that the selection of higher modes 
requires non-linear deformations from some (unstressed) 
reference configuration, we speculate that the mode num- 
ber n will be reducedby this effect. However, this assumes 
the increase in crust thickness to be uniform. In reality, 
it is expected that the rate of deposition of new material 
to the shell will be higher for regions that are convex rel- 
ative to the solution (that is, that 'bulge' into it) than for 
concave regions, since convex regions will be exposed to 
a greater volume of solution and thus solute. This sug- 
gests a two-way coupling between deformation and crust 
thickening, a potentially rich problem for which further 
study would be desirable. 
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When comparing our work to the experiments of 
Pauchard et al. @, IE 0; Gorand et al. p[ and Kajiya 
at al. @ , a difference in notation should be observed: in 
these works the term 'instability' is used to refer to the 
formation of any mode n > 1, whereas we reserve the 
term for snap-buckling events (also symmetry-breaking 
deformations). From our findings, modes with n > 1 
(such as mexican hats) form when the characteristic de- 
formation is of the order of the thickness of the shell, 
and thus will be delayed with respect to the formation 
of the crust itself. Nonetheless the sequence of deformed 
shapes with increasing 9$ observed in p| appears to be 
consistent with our work, although without any theory 
for the formation of the crust, it is difficult to provide any 
close correspondence with experiments. For instance, no 



n > 3 mode has been observed in experiments; this may 
be due to the emergence of asymmetric modes (which 
for uniform shells are known to occur by A ~ 8 [l5j . 
before the emergence of even the n — 3 'crater' mode), 
due to dynamic effects, or more simply that the required 
shell parameters are outside of the experimental param- 
eter window. 
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